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A modified phase field crystal model in which the free energy may be min- 
imised by an order parameter profile having isolated bumps is investigated. 
The phase diagram is calculated in one and two dimensions and we locate the 
regions where modulated and uniform phases are formed and also regions where 
localised states are formed. We investigate the effectiveness of the phase field 
crystal model for describing fluids and crystals with defects. We further con- 
sider a two component model and elucidate how the structure transforms from 
hexagonal crystalline ordering to square ordering as the concentration changes. 
Our conclusion contains a discussion of possible interpretations of the order 
parameter field. 



I. INTRODUCTION 



Modelling materials at the atomic scale is a task which, for example, may be performed 
using Molecular Dynamics simulations. This involves solving coupled equations of motion 
to calculate the position of each particle at every time step. The resulting calculations can 
be very computationally expensive, especially when one seeks to consider phenomena which 
involve a large number of particles. Only short atomic time scales can be feasibly accessed 
with this or other such approaches. However, there are some instances where it is important 
to consider materials on the atomic length scale for much longer diffusive time scales, e.g., 
when investigating freezing or glass transitions. One approach to such problems that may 
be adopted is to develop a phase field model capable of describing the structure of materials 
on the scale of the individual particles. In contrast to traditional phase field models, the 
recently developed phase field crystal (PFC) models are capable of just such a description 



2 



and are now widely used in the literature to model crystalline structures [TH5]. The PFC 
model consists of a Swift-Hohenberg-like equation j6J, but with conserved dynamics rather 
than the non-conserved dynamics of the regular Swift-Hohenberg equation. Similar models 
with conserved dynamics arise in quite different contexts as well [71 [8]. The regular PFC 
model is governed by the following equation: 

g0(x,t) _ 

dt 50(x,t)' ^ ^ 

where the free energy functional 



F[<P] = J f{<p), (2) 

with 

m = ^[r + {<l' + Vy]<P+^, (3) 

where a is the mobility coefficient, r is the undercooling parameter that decreases with 
decreasing temperature, g is a constant which determines the typical microscopic length 
scale in the system and (/)(x, t) is the order parameter. For certain parameter values this 
free energy functional is minimised by an order parameter profile consisting of a periodic 
array of bumps which somewhat resembles the density distribution of particles in a crys- 
talline material. This interpretation is bolstered by the fact that it has been shown that 
the PFC model (Eqs. Q and (g) may be derived from the density functional theory of 
freezing [9] and the dynamical density functional theory for colloidal particles [lOl [11] with 
certain approximations. The free energy is minimised by either periodic structures or by a 
homogeneous fiat profile, depending on the values of q, r and (f) = j dyi </)(x, t), where L'^ 
is the size of the system. In two dimensions (d = 2), one observes a homogeneous phase, two 
hexagonal phases (hexagonally ordered bumps/holes) and a stripe phase |TI43t[T2]. The liter- 
ature largely focuses on the region of the 2d phase diagram which contains the hexagonally 
arranged bumps and their transition to the homogeneous state [IHSl [9l [TOl [12] . The uniform 
profile 0(x, t) = (f) represents the order parameter in a uniform liquid and the hexagonal 
phase is treated as a crystal. The model is then used to consider a number of problems 
including melting and freezing [Sj [9l [TD] and grain boundary effects [HIS], [1]. 
In the 'standard' PFC model (Eqs. Q and (|2])) the hexagonally arranged bumps are con- 
sidered to be particles/colloids in a crystalline structure. The interpretation of the striped 
and hexagonally ordered hole structures is unclear and as such these phases are commonly 
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ignored. The conjecture that the ordered bumps represent crystaUine particle structures can 
be extended by including a 'vacancy term' in the free energy [T3| [T3] which strongly breaks 
the hole-bump (0 — )■ —0) symmetry of Eq. ([s]): 



+ /.ac(0) 



(4) 



Using this free energy Q, it is possible to obtain structures which contain a mixture of 
bumps and vacant areas (areas where the order parameter is approximately uniform around 
the value ~ 0), which in the interpretation of Ref. [I3] resemble snapshots of fluid con- 
figurations or crystalline structures with defects. We will return to the issue of the precise 
interpretation of the nature of the order parameter field in the conclusion. In this paper 
we investigate the thermodynamics and the structures formed in this augmented conserved 
Swift-Hohenberg model, or 'vacancy phase field crystal' (VPFC) model, and also in a two 
component generalisation of this model. The vacancy term takes the following form [T3l ITi] : 

UM) = Hct>Wcl>\-<t>), (5) 

where if is a constant. We use the value H = 1500, as in Refs. [131 El- This acts as a 
piecewise function which is zero for positive values of and takes an increasingly large value 
when < 0. Hence, this term penalises negative values of 0. This leads to the VPFC model 
forming periodic structures which are somewhat different from those of the regular PFC In 
addition, the VPFC model has a large region of parameter space at small 0, where spatially 
localised structures form. The time evolution of the order parameter is governed by the 
conserved dynamics used in the standard PFC model Q. 

We begin in Sec. |TT] by considering the phase behaviour of the model, investigating the 
transition between periodic and localised states. We focus on understanding the bifurca- 
tion diagrams connecting the various uniform, periodic and localised states exhibited by the 
model. We then go on to consider how individual localised states or particles interact with 



one another. In Sec. Ill we extend the model to consider a two component system, and we 
determine how the particles in the binary mixture interact with one another. We find a tran- 
sition between hexagonal and square ordering of the particles as the concentration changes. 



Our conclusions follow in Sec. |IV[ and include a discussion of the proper interpretation of 
the order parameter field 0. 



4 



II. ONE-COMPONENT SYSTEM 

A. Linear stability of a homogeneous profile 

We begin by considering the phase behaviour of the VPFC model (Eqs. ([T| and (|4])). We 
calculate the limit of linear stability for a homogeneous flat state using a linear stability 
analysis. In the context of colloidal suspensions exhibiting microphase separation and fluids 
of charged particles, this limit of linear stability is referred to as a 'A-line' |15HT8]. Since 
fvac is non differentiable at = we treat it in a piecewise manner, by treating the two 
cases (j) > and (p < separately (in fact, if 0(x) takes the form of Eq. ^ and (p > |^|, 
then fyac = everywhere and the thermodynamics of the VPFC model reduces to that of 
the regular PFC model). We assume that the order parameter takes the form of a flat 
profile plus an additional small amplitude harmonic modulation: 

= + 5/. = + ee^'^-'e'^*, (6) 

where (p is the average value of the order parameter and the amplitude |^| ^ 1. Substituting 
this expression into the functional derivative of the free energy Q we obtain: 

_ = (r + g^)0 + 3Hm\ - 0) + 0' + [{k' - q'f + A]&P + Oi&p'), (7) 

where 

A = r + 6i7(|0|-0) + 30l (8) 

Inserting this expression for the functional derivative ([T]) into the dynamical equation ([T]) 
and then linearising we arrive at the following dispersion relation: 



f3 = -PaHe - + A]. (9) 

When the growth rate P{k) > 0, any small amplitude modulation with wave number k = |k| 
will grow over time. There is a local maximum in /3 (which becomes the global maximum 
when the uniform state is unstable) at the wave number: 

= ^i/6g2 + 3v^4_3^_ (10) 

Thus, if one takes an initially almost flat profile 0(x, t = 0) = + A:'(x), where A'(x) is 
composed of a sum of a large number of small-amplitude harmonic modulations [cf. Eq. (|6])] 
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FIG. 1: Dispersion relation curves for the VPFC model (Eqs. ([TJ and Q) when g = 1. Four cases are 
shown, with (i) /3(fcm) > (red solid line), (ii) i3{km) — (blue dashed line), (hi) P{km) < but fc™ > 
(green dotted line) and (iv) /3 < and km = (magenta dash-dotted line). 

with different wave numbers k (in practice X{x.) is generated by adding a small random 
number to the discretised initial profile), then as the system evolves in time 0(x, t) will 
develop spatial modulation on the length scale since this scale corresponds to the maxi- 
mum growth rate Pm = P{km)- This length scale has an inverse dependence on the value of 
q, i.e., increasing the value of q reduces the length scale of the structures which are formed. 
The limit of linear stability is defined as the locus of points in parameter space where the 
maximum in the dispersion relation ^ is at zero, i.e., (3m = 0. The conditions /3 = = 0, 
subject to the requirement that ^ yield A = 0, /c^ = ±q. Thus A in Eq. ([s]) can be 
considered as a measure of stability: when A < the system is linearly unstable and when 
A > the system is linearly stable. The magnitude of A indicates how 'far' we are from the 
limit of stability. Figure [l] shows the dispersion relations (3{k) for various values of A. In 
accordance with Eq. ( [I0| ) the maximum aX km ^ q disappears when A > y; in this case only 
the maximum at A; = remains. It is important to note that these results are identical to 
those of the regular PFC model (Eqs. ([T]) and ([2])) for positive values of the order parameter 
> 0. 
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B. One-dimensional model 

In order to develop a better understanding of the effect of the 'vacancy term' (|5| we initially 
consider the phase diagram for the system in one spatial dimension. The regular PFC model 
(Eqs. ([!]) and (|2|) in one dimension exhibits two distinct phases pj: a non- uniform state 
in which the order parameter profile resembles a sinusoid and a uniform state in which the 
order parameter is a constant. The phase diagram of the regular PFC model is symmetric 
around = owing to the symmetry of the free energy (|2| with respect to — )■ —0. This 
is no longer the case when the vacancy term ([s]) is added. 

The phase diagram for the ID VPFC model is shown in Fig. |2](a) and is very different from 
that of the regular PFC [2J. As with the regular PFC model, modulated profiles are present 
below the limit of linear stability A = (blue dashed line) provided < \J^/2q^ . However, 
with the added vacancy term ([s]) the lower limit for the presence of the modulated phase is 
at > {H ^ 1). The tricritical point with > (red dot) familiar from the PFC model 
remains. Above this point the phase transition between the periodic and homogeneous 
phases is of second order. Below this point a periodic phase with = 0p coexists with a 
homogeneous phase with = 0^ and the phase transition between these phases is of first 
order. Figure [2]^a) shows the coexisting phases using fixed temperature (horizontal) tie-lines 
connecting 0p and (ph (solid red lines). 

We can calculate the location of the tricritical point as follows: Since the wavenumber near 
A = is ~ Q' we assume that the order parameter profile takes the form 

(j) = (pp + A cos qx + B cos 2qx + . . . (11) 

and compute the free energy F. When 0p > the vacancy term drops out and we obtain 
the following expression for the free energy per unit length fp = F/L of the periodically 
modulated phase: 

(12) 



We refer to the Ansatz (11) as the two-mode approximation. This approximation is reliable 
around and above the tricritical point since the amplitude of the modulation in is small 
when |A| ^ 1. Moreover the two-mode approximation appears to be exact at the tricritical 
point, since the location of the tricritical point is unaffected by the inclusion of cosSgx 
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FIG. 2: The phase diagram for the ID VPFC model (Eqs. ([I]) and Q) is displayed in (a) for the case 
q ~ I. The red solid lines are the coexistence curves between the periodic and uniform phases; the red 
circle is the tricritical point. The blue dashed line is the locus A = 0, which is the limit of linear stability 
for uniform profiles. The green dot-dashed lines are a guide showing the parameter space where local and 
periodic structures are formed, (b) - (g) show examples of order parameter profiles from numerical simu- 
lations corresponding to (local) minima of the free energy, for the values of (f> and r indicated in (a). The 
parameter values are: g = 1, r = —0.9, a — 1 and (b) (f> — 0.01, (c) (f> — 0.05, (d) 4> = 0.1, (e) (j) — 0.175, 
(f) ^ = 0.3 and (g) ^ = 0.5. 



and other higher order modes. In contrast, the mode Bcos2qx must be retained in order 
to obtain the correct value of the amphtude A in the vicinity of the tricritical point (see 
below). 
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To demonstrate this we minimise fp in Eq. (12) with respect to the amphtudes A and B, 
obtaining the following two conditions 



r + 30$ + 30p5 + ^A^ + ^B^ = 0, 



(r + 9q')B + 3^lB + l^pA' + ^A^B + 



0. 



(13) 
(14) 



Solving these for the amphtudes A and B, and substituting into Eq. (12), we obtain an 



approximation for the free energy density of the periodic phase fp. Linearising Eq. (14) in 
B, we find that 

B = -^ + 0{ApA\A') (15) 
where Ap = r + 30^ and hence, from Eq. (13), that 



A 



(16) 



The free energy density of the homogeneous phase fh, having (f){x) = (ph, is obtained simply 
setting A = B = in Eq. (12) to obtain fh = h{r + q'^)4>1 + The chemical potential in 



the homogeneous phase is = dfh/dcph, and in the periodic phase fip = dfp/d(j)p. 

To calculate the location of the tricritical point we recall that at coexistence between the 

periodic state and the homegeneous state we must have ^ip = fih- We write the average value 

of in the periodic state ipp = 4>h + C, where C is the difference between the average value 

of the order parameter in the two coexisting phases, implying that the coexistence condition 

is 



f^p{4>h + C) - fih{4>h) = 0, 



(17) 



or equivalently: 



3^ 



(r + q^)C + 3^lC + -^hA^ + 0{A\ CA\ C^) 



0. (18) 

Since the amplitude A of the modulated phase at coexistence is small when \Ah\ ^ 1, where 



Ah = r + 30^, Eqs. (14) and (18) yield, for \Ah\ <^ 1, the expressions 



B 



(phA^ 



+ 0{AhA',A^' 



C 



rA^ 



2g4 



+ 0{AhA\A'). 



(19) 



Equation (13) then yields 



A = 2 



A,/ _380| 
3 V 3g4 



-1/2 



+ 0(A,,). 



(20) 
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Taking the limit C — )■ now takes us to the tricritical point. At the tricritical point 
Ah = Ap = and the chemical potentials iip{(t>h) and iihi.'Ph) are identical. Thus from Eq. 
(20) we see that the tricritical point occurs at = ^/S/SSg^ ~ 0.281g^, r = — (9/38)g'^ ~ 
— 0.237g^. These coordinates agree precisely with the result obtained from the common 
tangent construction between the free energy of the periodic phase and the free energy of the 
homogeneous phase to determine phase coexistence, and also with our numerical simulations 
of the VPFC model performed with q = 1. Thus when g = 1 the phase transition is of second 
order for r > —9/38 and of first order for r < —9/38. As already mentioned this result is 



exact in the sense that it is unchanged if more modes are included in the Ansatz (11) 
and improves on the prediction r = —1/4 obtained for the PFC model using a one-mode 
approximation [2j. Note that the vacancy term ([s]) does not affect the transition because 
is positive everywhere. 

As discussed above the transition between the periodic and the uniform phases is of first order 
below the tricritical point. As r decreases, this region of the phase diagram is increasingly 
affected by the vacancy term ([5]), as the amplitude of the structures becomes large enough 
to reach negative values. We observe that including the vacancy term ([s]) decreases the 
distance between the coexistence curves. This is because the vacancy term increases the free 
energy of the profiles in the periodic phase, which decreases the difference between the free 
energy of the periodic structures and the homogeneous state and hence a common tangent 
construction between the two yields values which are closer to the linear stability line A = 0. 
We calculate the coexistence values below the tricritical point by numerically solving for the 
order parameter profiles, because the two-mode approximation becomes inaccurate when 
r < —0.3, where the order parameter profile develops regions where < and so the 
vacancy term makes a contribution to the free energy. Note that for the regular PFC model 
(i.e., when H = 0) the two mode approximation for the free energy works very well, agreeing 
to two significant figures or more with the exact free energy for r > —0.9. 
To calculate the coexisting phases numerically we select a value of r and determine the order 
parameter profile along this line for different values of 0. The free energy at each of these 
points is minimised with respect to the domain size, which effectively gives us the minimum 
free energy for the infinite system. A polynomial is then fitted to these values to produce a 
continuous curve which gives the free energy of the periodic phase for the chosen value of 
r. We then make the common tangent construction between the free energy of the periodic 
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phase and the uniform phase to calculate the two coexisting values at the chosen value of 
r. This process is then repeated for different values of r. The resulting coexistence curves 
are plotted as the solid red lines in Fig. ^a). 

The periodic structures which are formed by the VPFC model [Fig. |2](f)] are qualitatively 
very similar to the structures which can be found in the regular PFC model. However, the 
amplitude of the modulations is restricted by the large penalty in the free energy accumulated 
when (f) < 0. Inside the coexistence region between the periodic and uniform states, we 
observe interesting structures where the amplitude of the peaks does not remain constant 
and a second length scale is visible in the structures [Fig. |2]^g)]. This is also an effect 
which is present in the regular PFC model and will be discussed in detail in future work. 
What is most intriguing, and is perhaps the most appealing aspect of the VPFC model, 
is the appearance of localised states for small positive values of (f) when the magnitude 
of r is sufficiently large (r < —0.6). We obtain order parameter profiles by numerically 
integrating forward in time Eqs. ([T]) and Q until a stationary solution is reached, starting 
from the initial profile 0(x, t = 0) = + X(x.), where X is a small amplitude random noise 
profile with zero mean. A rich variety of different patterns is observed, including periodic 
structures mixed with almost flat regions [Fig. ^d)] and individual isolated peaks [Figs. [2]^b) 
and (c)]. In Fig.|2ja) the green dot-dashed curves indicate the boundary of the region where 
one observes regular periodic structures and where the localised structures are formed. Note 
that these are guidelines only and are not thermodynamic coexistence curves. The lower-left 
dot-dash curve roughly denotes the linear stability limit of the regular periodic structures, 
such as that in Fig. ^f). This is determined numerically. We begin with a periodic proflle 
and reduce the value of gradually, minimising the free energy at each step, while keeping 
r constant. The limit point is then deflned as the value of where the periodic proflle 
becomes linearly unstable and a vacancy is introduced. In a similar way, we determine the 
upper-right dot-dash line, which is the limit of linear stability of the structures with defects. 
This is found by starting with a proflle containing a single vacancy and increasing until the 
vacancy disappears. These two points are calculated for different values of r and then a best 
flt to this data is shown in Fig. [2]^a). There is some hysteresis in the region between these 
two curves, with the type of proflle produced depending heavily upon the initial conditions. 
Within the localised state region of the phase diagram it is possible to obtain order parameter 
proflles with a varying number of peaks for a given system of length L. Keeping r < —0.6 
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constant and varying allows us to control the number (density) of bumps as shown in 
Figs. [2](b)-(e). Beginning with = we find isolated peaks in large vacant areas (where 
is approximately uniform with 0^0). As is increased the number of peaks increases 
until we return to the familiar regular periodic structures. The assumption of Ref. [I3] is 
that unlike in the regular PFC, where the uniform phase is associated with the liquid and 
the modulated phase with the crystal, in the VPFC model one may associate each bump in 
0(x) as corresponding to a particle and so the model can describe fiuids [Figs, and (d)], 
crystals with vacancies and defects [Fig. ^e)] and regular crystals [Fig. |2](f)]. 
The findings presented in Fig. [2] indicate the existence of a hysteretic transition between 
periodic and localised states, and are a consequence of homoclinic snaking [T9H22] in the 
present system. In the standard homoclinic scenario such localised states are present within 
a part of the coexistence region called the pinning region. The localised states in the lower 
left part of the parameter plane (0, r) in Fig. [2]^a) correspond to the global energy minimum 
or to other deep but local energy minima. Families of such steady state solutions can be 
obtained for the VPFC model that we study here [Eq. ([T]) with Eqs. (|3]), ^ and ^] by 
employing the path continuation techniques bundled in the package AUTO07p p3]. As an 
example, in Figs.|3]and|4]we show the characteristics of localised solutions along cuts through 
the plane (0, r). In particular. Figs, [s] and |4] give results for changing r (at constant = 0.1) 
and changing (at constant r = —0.9), respectively. All solutions are characterised by their 



norm ||50|| = y (1/L) {(p{x) — (pydx, chemical potential /i = 6F/6(j), mean free energy 
density difference (-F[0(a;)] — Fq)/L, where Fq = F[0] and mean grand potential density 
u = F[(f){x)]/L — (f)fi, and satisfy periodic boundary conditions on the domain < a; < L. 
It turns out that there exist three types of localised steady states: (i) the heavy solid black 
line consists of a; — > —x symmetric localized states that have a maximum at the centre, 
i.e., the overall number of bumps within the structure is odd. (ii) The dashed red line also 
represents x — —x symmetric localized states but this time with a hole (minimum) at the 
centre, (iii) The localized solutions of the third type are not symmetric under x — > —x 
and are called "asymmetric states". These reside on branches that connect (via pitchfork 
bifurcations) the two branches of symmetric localized states. These branches are included 
in the bifurcation diagrams as dotted green lines. 

Examples of order parameter profiles of types (i)-(iii) are presented in Fig. [sj corresponding 
to the various solution branches displayed in Fig. |4} This sequence of profiles expands upon 
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(c) - - - r - (d) 

FIG. 3: (color online) Bifurcation diagram showing localized solutions of the VPFC (an augmented con- 
served Swift-Hohenberg equation) [Eqs. ([l]) and with H = 1500, as a function of the parameter r, for 
the mean order parameter = 0.1 and a fixed domain size of i = 100. The various solution profiles are 
characterised by their (a) norm, (b) chemical potential /i, (c) mean free energy density {F — Fq)/L, 
and (d) mean grand potential density lu = F/L — (f)fi. The heavy black dash-dotted line corresponds to 
the homogeneous solution (j){x) — cj). Periodic solutions with n = 15 bumps are shown as a thin blue 
dashed line, whereas the nearby thin black dotted lines represent the n — 14 and n = 16 solutions as in- 
dicated in the plot. The heavy solid black and dashed red lines that bifurcate from the n = 15 periodic 
solution represent symmetric localized states with a maximum (odd states) and a minimum (even states) 
at the center, respectively. The green dotted lines that connect the two branches of symmetric localized 
states correspond to asymmetric localized states. Together the branches of localized states form a slanted 
snakes-and-ladders structure. 



the few examples shown in Figs. [2|^b)-(g). Recall, however, that the results in Figs. |2]^b)-(g) 
are obtained starting from an order parameter profile with a small amplitude random noise 
and so they do not always exactly agree with the steady states at the same resulting 
from the path continuation. The norm, chemical potential /i, mean free energy density 
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FIG. 4: (color online) Bifurcation diagram showing localized solutions of the VPFC as a function of the 
mean order parameter 0, for r — —0.9 and a fixed domain size of L = 100. The various solution profiles 
are characterised by their (a) norm, (b) chemical potential /i, (c) mean free energy density {F — Fq)/L, 
and (d) mean grand potential density uj = F/L — (f>^. The line styles are as in Fig. [sj Here, however, 
the heavy solid black and dashed red lines bifurcate at large (j) from the n = 15 and n — 14 periodic 
solutions, respectively. Typical profiles for all the branches of localized states are given in Fig. [5] The 
vertical dotted lines in (a) correspond to values of for results in Fig. [s] The blue dots correspond to the 
five time simulation profiles shown in Fig. j2|^b) - (f). 

{F — Fq)/L and the mean grand potential F/L — 0/i have been calculated for the profiles 
obtained from time simulations [Figs. [2]^b)-(f)] and are plotted as blue dots in Fig. |4j A 
close inspection reveals that the energy of the time simulation results is often slightly higher 
than that from the continuation results, indicating that in these cases the time simulation 
converges to a local and not the global energy minimum. This is to be expected as the 
solutions shown in the bifurcation diagrams are only the 'tip of the iceberg'. For instance, 
there exist many more solutions, where not all the inner distances between the bumps are 
identical. This is related to the fact that individual bumps have oscillatory tails and the 
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FIG. 5: A selection of steady state profiles (t>{x) for r — —0.9 at = 0.01, 0.05, 0.1, and 0.2. From 
top left to bottom right we show first nine type (i) solutions, i.e., symmetric localised states with an odd 
number of maxima (in black), then eight type (ii) solutions, i.e., symmetric localised states with an even 
number of maxima (in red), followed by six type (iii) solutions, i.e., asymmetric localised states (in green). 
The final image is the n = 15 periodic solution at (/> = 0.3 (in blue). The number in each panel indicates 
the corresponding value of <j}. The solutions from the symmetric branches are shown in the sequence that 
follows the respective branch in Fig. |4]ja), starting from the left. The asymmetric states for identical (j> are 
shown in the order of decreasing norm. 



'locking of these tails' allows for different equilibrium distances [21]. The solutions presented 
in Figs. |3] and |4] represent the solution having the lowest energy in the respective class. 
However, the energy differences between these and the 'less symmetric' solutions are often 
tiny. Thus, it is not surprising that time simulations starting from random initial profiles 
often converge to solutions with greater disorder and energies than those shown in Fig. |4} 
For instance, the solution in Fig. ^^c) at = 0.1 is a nine-bump solution similar to the 
odd symmetric localised states shown in the first two panels of the second row of Fig. |5} 
The amplitudes agree well and although the arrangements of the nine bumps are different, 
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the free energy and norm still agree to < 1%. However, at large average order parameter 
values the time simulation results can converge to metastable states with energies quite 
different from the minimum energy states for domains of this size L = 100. For example, 
the periodic solution obtained from the time simulation (shown in Fig. ^i)) when = 0.3 
has eighteen bumps. However, from Fig. |4](c) we observe that the energetic minimum is 
obtained by a periodic profile with fifteen bumps, as shown by the steady state solution 
in Fig. |5| The convergence to a different number of bumps in the time simulation may be 
caused by discretisation effects or by the initial noise profile used. As one would expect, the 
free energy associated with the eighteen bump periodic structure is significantly larger than 
the fifteen bumped profile. 

In Fig. [3] (0 = 0.1) the localised states bifurcate subcritically from the periodic solution 
branch (that itself emerges from the trivial homogeneous solution that is displayed as the 
heavy black dotted line). Therefore, one expects hysteretic behaviour as encountered in the 
time simulations. A magnification (not shown) allows us to determine the threshold values 
for the hysteretic transition. When decreasing r in the region where periodic solutions are 
always found, one first passes Vsn = —0.685 where the last 2 branches of localised solutions 
annihilate in a saddle-node bifurcation [Fig. |3|^a)]. Slightly below r^n, both the periodic 
solution and the localised state with a single bump are local energetic minima. Although 
the periodic solution represents the global minimum, particular time simulations sometimes 
converge to the localised state. The differences in energy between the two is < 1% in the 
case of Fig. [3j When r is further decreased below r^^ = —0.700 the energy of the even 
symmetric states becomes smaller than the one of the n = 16 periodic solution, that is 
however still linearly stable. The situation changes at Vc = —0.749 where both symmetric 
localised branches bifurcate from the n = 16 branch, i.e., below Vc the latter is linearly 
unstable. Furthermore, below the energy of all localised states rapidly becomes much 
smaller than the energy of all periodic states [Fig. [3|^c)]. The hysteresis range displayed in 
Fig. [2]provides a good approximation for the region between Vc and rsn. This region becomes 
larger as is increased. 

The situation is very similar when is changed for fixed r (Fig. |4|. The resulting hysteresis 
range is between = 0.150 and 0.239 for symmetric localised states with an odd number 
of maxima and between = 0.202 and 0.265 for symmetric states with an even number 
of maxima. Overall, one should therefore expect a wide hysteresis region roughly between 
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= 0.15 and = 0.25. The hysteresis range obtained from the time simulations (indicated 
in Fig. [2| is roughly 0.19 < (p < 0.22. This is narrower than the range deduced from the 
path continuation analysis of the localised steady states, but lies right in the middle of it. 
Before we move on to discuss the two-dimensional case, we should comment on how our 
results fit into the wider context of research on localised states. Much research on localised 
states focuses on the non-conserved Swift-Hohenberg equation [T9H2T]. There, such states 
can only exist if the primary bifurcation of periodic states from the homogeneous base state 
is subcritical. The localised states exist in a sub-range of the existence range of the periodic 
states bounded on either side by the saddle-node bifurcations of the branches of symmetric 
localised states. In the non-conserved Swift-Hohenberg equation these accumulate expo- 
nentially rapidly towards the parameter values corresponding to first and last tangencies 
between the unstable manifold of the homogeneous state in space and the stable manifold 
of the periodic state. These tangencies define the pinning region containing the different 
localised structures. In contrast, in the presence of a conserved quantity localised states 
may exist outside the existence region of periodic states, may occur even in the supercritical 
case and the saddle-node bifurcations of the localised states are no longer aligned, i.e., one 
finds slanted snaking [25] . This is typically a finite size effect [^ . 

For the regular PFC (conserved Swift-Hohenberg equation) [Eq. ([T| with Eqs. ^ and ([s])], 
localised states are briefly mentioned in Ref. [8]. However, no systematic results along the 
lines of those presented in Refs. [IH] for non-conserved or [25] for conserved order parameter 
fields are available. The model used here is a special case because it includes the non- 
analytic vacancy term ([5]). However, a similar bifurcation structure is found for the classical 
conserved Swift-Hohenberg equation i.e., the regular PFC model [2j]. 

C. Two dimensional model 

We now move on to consider how the VPFC model behaves in two dimensions. As with 
the regular PFC model [2], when we expand into two dimensions we observe stripes [see 
Fig. |6](b)] and hexagonally ordered bumps or holes [see Fig. ^c)]. In Fig. [6|^a) we display 
the phase diagram of the VPFC model in two dimensions and typical time simulation results 
from the striped [6|b) and hole [6|c) phases, calculated on a regular grid with grid spacing 
dx = 0.5. Square ordering of bumps or holes does not appear in the phase diagram because 
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FIG. 6: The phase diagram for the 2D VPFC model (Eqs. ([T]) and Q) is displayed in (a) for the case 
q — I. The red solid lines are the coexistence curves between the various phases. The blue dashed line 
is the limit of linear stability for uniform profiles A = 0. The green dash-dotted lines indicate the re- 
gion where localised and hexagonally ordered bump structures coexist. Simulations of (b) stripes and 
(c) hexagonally ordered holes are also shown. The parameter values for these simulations are: q = 1, 
r = -0.9, a = 1 and (b) = 0.4 and (c) 4> = 0.53. 



these structures always have a higher free energy. However, this can be changed through 
appropriate alterations to the free energy [5] . Square ordering can also occur when extending 



to a two-component mixture (cf. Sec. Ill C below). Using the same method as outlined above, 



we calculate the regions of the phase diagram where there is coexistence between hexagonally 
ordered holes and the uniform distribution, between holes and stripes and between stripes 
and hexagonally ordered bumps. The vacancy term ([s]) shifts the modulated phases into 
the positive > plane. The section of the phase diagram where holes are observed is 
much smaller when compared to the regular PFC model and now extends beyond the limit 
of linear stability of the flat state (at A = 0). This means that for certain values of (where 
< A ^ 1), hexagonally arranged holes are energetically favourable but are only observed 
in time simulations for certain initial conditions - i.e., when starting with an order parameter 
profile 0(x, t = 0) which already has modulations which are sufficiently large in amplitude. 
As r is decreased (i.e., for larger |r|) it becomes increasingly difficult to obtain structures 
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with holes up to and inside of the coexistence region between the hole and the uniform 
phases. This is a direct consequence of the limit of linear stability occurring in the middle 
of the hole phase. Therefore, the accuracy of results for the coexistence region between 
the hole and uniform phases decreases as |r| becomes larger. The stripe phase occurs in 
between the two hexagonal phases. In the simulation order parameter profiles displayed in 
Fig. ^h) and (c) we observe various defects and in (c) 'grain' boundaries between regions 
with different orientations, which depend on the initial conditions (our initial profile was a 
flat state with additional small amplitude white noise). The true minimum profile for case 
(b) is a series of parallel stripes which are identical to the periodic profiles in the ID system 
(shown in Fig. [2|f)). 

The most important portion of the phase diagram from the materials modelling point of 
view, is the bump phase because the basic assumption is that each bump represents a 
particle. When r > —0.4 or when has a value close to that in the coexistence region 
between bumps and stripes, we observe hexagonally arranged bumps, similar to those in the 
regular PFC model. However, in a similar manner to the ID system, we observe localised 
structures at small values of when r < —0.4. In the phase diagram [6|a) the green dot- 
dashed lines are numerically obtained estimates for the location in the phase diagram of 
the limits of linear stability of the uniform periodic states (lower curve) and the localised 
(vacancy) states (upper curve). They are determined in the same manner as discussed above 
for the one dimensional system for a square system of side length L = 25. It is important 
to note that the parameter range where localised bumps coexist with regular periodically 
ordered bumps is much broader for the 2D system, implying a large amount of hysteresis. 
We now focus our discussion on the portion of the phase diagram where isolated bumps form. 
As can be seen in Fig. [7} these profiles resemble particle configurations in gases, liquids and 
crystalline solids and so the VPFC may be a valuable model for describing materials on 
the microscale [13] . This region of the phase diagram is full of complexity and many varied 
structures may be observed. However, here we forego a full systematic study of this large 
region in parameter space and limit ourselves to showing representative results obtained 
for a single value of the undercooling parameter r = —0.9 for which there is a fairly large 
range in with isolated bumps. We set the initial order parameter profile to be a uniform 
state with a small amplitude noise 0(x, t = 0) = + X{Y — 0.5), where F is a random real 
number uniformly distributed between and 1 and A = 10^^ is the amplitude of the noise. 
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FIG. 7: (a)-(c) Typical steady state order parameter profiles obtained in time simulations for increasing 
(j). (d) - (f) show the corresponding radial distribution function g{x) calculated from multiple simulations. 
The parameter values are: a — 1, q = 1, r = —0.9 and in (a), (d) (j) — 0.01, in (b), (e) (f> — 0.1 and in (c), 
(f) 4> = 0.24. 



We consider three cases; i) (f) = 0.01 and ii) (j) = 0.1 where a disordered arrangement of 
locahsed bumps forms and iii) = 0.24 which is in the region where bumps are hexagonally 
ordered. We average over many simulations to calculate the two point correlation function 
for each of these cases. This is done by locating all the maxima in the equilibrium profile 
(/)(x), for a given initial realisation of the noise; i.e., we locate the position of all the bumps. 
From these sets of coordinates we calculate the radial distribution function g{x) in the usual 
way [28]. We display a simulation result for case i) in Fig. [Tj^a) and the corresponding 
radial distribution function g{x) in Fig. [7](d). We find that there is almost no correlation 
between the bumps in this circumstance except for the core repulsion and a very small 
peak at a; ~ 16, indicating that there is a weak attraction between the bumps. Therefore, 
simulations with these parameter values appear to qualitatively describe gas-like formations 
of particles/colloids. In Fig. ^h) and (e) we plot a typical order parameter profile and 
the corresponding g{x) for case ii). We observe a large increase in the number of bumps 
as compared to the previous case. The radial distribution function shows that we have 
strong short range ordering, but without any long range order. This is very reminiscent of 
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the ordering in liquids. There is a very sharp peak in g{x) at around x = 7.5 (which is 
approximately the diameter of the bumps) and a smaller peak around a; = 15. A similar 
example is also given in Ref . [13] . If we further increase the value of we eventually find the 
more familiar hexagonally structured array of bumps which is reminiscent of the ordering in 
simple crystalline solids. In Fig. [7](c) we display an example of the order parameter profile 
for case iii) and in Fig. [7]^f) we show the corresponding g{x). For this case we observe that 
g{x) is highly structured indicating the system has very strong short range correlations with 
a significant degree of long range ordering. We observe the split second and third peaks, 
which is a classic sign of crystalline order. These results indicate that the VPFC model 
may be used to model crystalline structures, much like the regular PFC model. The major 
difference between the two models is the existence of the fluid-like configuration of bumps 
observable in the VPFC model. In contrast, the fluid phase in the PFC model corresponds 
to the homogeneous state. 

The variation in the size and shape of the bumps that are formed is fairly small. In Fig. [sj^a) 
we display a selection of results for the order parameter profile through the centre of the 
bumps for the case when r = —0.9 and = 0.01. We determine the shape of the bumps by 
plotting the value of the order parameter (p against the distance from the peak of each bump 
(as shown by the data points). We can then fit functions which take the following form: 



e{x) = Poe-^'"' "'^'^ "^'"^ cos(/34x) + /35. (21) 

We fit this form to the data using a least squares method. The exponential part of 6{x) 
describes the decay of the modulation as the distance from the peak increases and the cosine 
function captures the oscillatory tail of the modulations which is an important factor in their 
interaction with other bumps |29l|30]. Figure [s]^ a) displays two cases; the (+) points and 
red solid line show the case q = 1 and the (x) points and blue dashed line show the case 
where q = 1.1. The size of the bump is reduced as we increase the value of q. This is because 
increasing the value of q increases the typical wave number which results in a smaller typical 
length scale (cf. Fig. [l] and Eq. (10)). 

The curves obtained from fitting the bump profile can be used to obtain an approximation 
for the effective pair potential V{x) between two isolated bumps, where x is the distance 
between the centres of the bumps. We take a uniform system with the value of equal to 
that in the uniform areas between bumps found in simulations for = 0.01, corresponding 
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FIG. 8: (a) Several sets of numerical results for the order parameter profile through the centre of a bump 
(+) for q = 1 and (x) for q = 1.1, together with fits to the data (solid red and blue dashed lines). These 
fits are then used to calculate the effective pair potential between two bumps. These pair potentials are 
displayed in (b). The inset displays a magnification of the tails of Vix). The parameters values are: r = 
-0.9 and 4> = 0.01. 



to the results in Fig. [7]^a). We then impose upon this the profiles for two bumps using the 
fitted curves shown in Fig. |8](a). We vary the distance between the superposed bumps and 
calculate the free energy of the system. We assume thereby that the two bumps retain their 
shape when they are close, despite the fact that in reality the bump shapes become distorted 
as bumps are pushed close together. 

In Fig. [8|^b) we display the results for g = 1 (red solid line) and g = 1.1 (blue dashed line). 
We observe that there is a shallow minimum in the potential at the distance x ~ 7.5 when 
q = 1 and at x ~ 7 when g = 1.1 (see inset of Fig. [Sj^b)). The minimum is at a smaller 
distance when q is larger because of the decreased diameter of the bumps - recall that q 
determines the size of the bumps. The resulting weak attraction between the bumps may 
also be inferred from the radial distribution function g{x) calculated for the low density case 
= 0.01 when q = 1 displayed in Fig. [Tj^d). We observe a second minimum in the potentials 
at X ~ 3.15 when q = 1 and at x ~ 2.3 when q = 1.1, where the former rather appears 
like a 'shoulder'. The order parameter profiles at the second minima resemble the elongated 
almost elliptical shapes which are observed in and around the coexistence region between 
bumps and stripes. See also Fig. ^c) where we also observe elliptical holes along some of 
the grain boundaries. 
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III. TWO COMPONENT SYSTEM 

We now extend the model to consider a binary mixture in perhaps the most simple way 
possible, by adding together free energies like Eq. Q for two order parameter fields 0a (x, t) 
and 0f,(x, t). We introduce a simple coupling term which allows the two components to 
interact with each other. This gives us the following expression for the free energy: 



(22) 



/(0a(x, t)) + fvac{(l)a{^, t)) + /(0fe(x, t)) + /^,ac(0b(x, t)) + r](f)a(j)b , 

where r] is the coupling coefficient and the functions / and f^ac are defined as before in 
Eqs. ^ and ([s]). The value of r is set equal for both components. However, we allow the 
value of q to be different for each species, so we now refer to these values as qa and qb, where 
the subscript denotes the corresponding component. Setting different values for q in the 
two components (i.e., ^ q^) results in an asymmetrical system in which the size of the 



bumps/modulations in (pa differs from that in 0^, as discussed further below in Sees. IIIB 



and IIIC Note that a different coupling term is used in Ref. [11]; a somewhat different 
two-component PFC model is presented in Ref. [31] . 

Here, just as for the one component model, we assume the dynamics of the system is governed 
by the following pair of equations (cf. Eq. ([T])): 

5F 



dt 



■'a 



at 

We also assume that the two mobility coefficients are equal: = = «. The two 
components are coupled purely by the term ri(pa<Ph in the free energy. When 77 > 0, this 
coupling term leads to a repulsion between the two species and so penalises structures which 
overlap or form on top of each other. The value of the parameter 77 determines the 'strength' 
of the coupling, and so the two component model reduces to two disconnected one component 
models in the limit 77 — )■ 0. 

A. Phase behaviour 

When the coupling coefficient is fairly large ri > 0.1, the coupling term has a significant 
impact on the phase behaviour of the model. In particular, the limit of linear stability and 
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the phase coexistence curves extend to much larger values of = 0a + 06 than for the one 
component model. We now determine the linear stability of a flat state in the model. We 
assume that the order parameter profiles of both components take the form: 

0a = 0a + 5/) = 0a + ee^'"e^*, 

0, = 0b + x^ = h + xie'^'^eP'. (24) 

where the amplitude |^| ^ 1 and the parameter x is the ratio between the amphtude of 
the modulations in the two components. The sign of x indicates whether instabilities are 
in-phase (x > 0) or anti-phase (x < 0) between the two coupled order parameter fields. 
From the magnitude of x we can deduce whether the instability is initiated from species a 
(Ixl ^ 1), species h {\x\ ^ 1), or a combination of both = 0(1)). We make a Taylor 
series expansion of the functional derivatives of the free energy with respect to the two order 
parameters 0a and 0^, to obtain: 

5F - ----3- 

7— = {r + ql)(j)a + ^H(j)a{\(j)a\- 4>a) + (t>a + ITPb 
d(Pa 

+ [{e-qlf + ^a + Xv\^ + Om, 

6F - - - - - 3 - 

7— = (r + g^)0fe + 3if06(|0fe| - 0;,) + 06 +?70a 

+ [x{k' - qlf + xAfe + 77] 5/. + 0{5f). (25) 
where A is defined as before (Eq. ([s])) and the subscript denotes the corresponding compo- 



nent. We substitute these expressions into the dynamical equations (23), yielding the matrix 
problem 



H = M , (26) 



where 

M = -k'^a 



{ql - k^f + Aa 

V {ql - + Afc 

We can now determine the dispersion relation P{k) by calculating the eigenvalues of M: 



,(,).l£M±^l£M!_|M|. (27) 

The resulting dispersion relation (]{k) is a double- valued function. However, since the growth 
rate along the + branch is always larger than that along the — branch, the limit of linear 
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stability can be determined from the + branch alone. If we assume that qa = Qb 
dispersion relation simplifies significantly, yielding: 



q, the 



2{e - + A„ + A, - v/(Aa- A5)2 + 4r/2 



(28) 



There is a local maximum of this expression which occurs at the positive wave number: 



kn 



24g2 + 6 4g^ - 6(A, + A^) + 6v/(A„ - As 



4^2 



(29) 



Substituting this wave number back into the dispersion relation (28), allows us to calculate 
the parameter values such that P{km) = (i.e., the limit of linear stability of a fiat state). 
We arrive at the following relation: 

A,A, = r]^. (30) 



When the system is linearly unstable it is possible for P{k = 0) to be a minimum or maximum 
(this transition occurs at A = — in the one component model). This is equivalent to the 
coefficient of k'^ changing from a positive value (minimum) to a negative value (maximum). 
The sign of the coefficient of k"^ is determined by the sign of the following quantity: 



Co 



d 9 d^g 

(902 502 



dg 



(31) 



where g{(j)a, 4>b) = f{,4'a) + fvaMa) + /(06) + fvac{4>b) +V4>a4>b- When C2 is negative/positive 
P{k = 0) is a minimum/maximum, this relation also holds for asymmetric systems where 
qa 7^ qb- In figure [9] we display typical dispersion relations when = g?, = 1, A^ = A;, = A 



and 7] = A. We show the case when i) the system is linearly unstable and C2 [Eq. (31)] is 



negative (red solid line), ii) the system is linearly unstable and Eq. (31) is positive (blue 



dashed line), iii) the system is at the limit of linear stability (i.e., Eq. (30) holds) (green 
dotted line) and iv) the system is linearly stable (magenta dash-dotted line). We observe 
that when qa = qb, the typical wave number k^ qa as we approach the limit of linear 
stability A^Af, — r/^ — ). 0. In the more general case with g^ 7^ g^ the dispersion relation may 
have two maxima at positive values of k neither of which occurs at g^ and qb- In this case 
the stability boundary is defined by the vanishing of growth rate I3{km) = of the larger of 
the two possible maxima of /3. 



From Eq. (30) it is clear that depending on the value of the coupling coefficient 77, the region 
of parameter space where the system is linearly unstable can be greatly larger than that for 
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FIG. 9: Dispersion relation curves for the two component VPFC model (Eqs. (22) and (23)), when = 
qb = 1, Aa = At, = A and rj — 4. Four cases are shown: (i) /3{km) > with ^{k = 0) a minimum (red 
solid line), (ii) f3{kra) > and f3{k = 0) a maximum (blue dashed line), (iii) j3{km) = and (iv) (i{km) < 
(magenta dash-dotted line). 



the one component system. For example, picking the value rj = 4 when r = —0.9 and setting 
the average value of both order parameters to be equal (pa = 4'b = 4>, we find that the limit of 
linear stability increases from = 0.548 (for the one component case) to (p = 1.278. As one 
would expect, this also increases the region of the phase diagram where modulated structures 
are formed. Our focus here is on the regions of parameter space where bumps are formed 
as this is the regime relevant to modelling crystalline solids. However, before proceeding 
to this, we make a brief survey of some of the structures which may be observed for larger 
values of (pa and (pb which lie outside of the bump phase. For the parameter values r = —0.9, 
1] = 4, qa = qb = 1 and (pa = (pb = (p show in Fig. 10 a sequence of order parameter 
profiles with increasing 0, for values of (p that lie above the region where bumps are observed 
(see later sections for a detailed analysis of the bump structures found in the two component 



model). In Fig. 10 we display scaled plots of order parameter profiles which are stationary 
states obtained from time simulations. We plot an order parameter defined as the normalised 
difference between the 0i(x) values of the two components A0(x) = (pa{'x)/(pa — (pb{'x.)/(pb 
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FIG. 10: The scaled order parameter A(/) for the two-component model, corresponding to minima of the 
free energy. The peaks in species a are shown in orange, peaks in species b are shown as blue and white 
areas show regions where (f>a ~ (t>b- The parameter values are: rj — 4, r = —0.9, Qa = Qb = 1, (pa — <Pb — 4>^ 
where (a) 4> = 0.25, (b) 4> = 0.3, (c) 1, (d) <^ = 1.15, (e) 4> = 1.2 and (f) 4> = 1.27. 



where 



j^max J^max 
V^a Vb 



0max 



(32) 



and where (p^^^ 



and 



are, respectively, the maximum and minimum values of 0i(x). Acj) 
is defined so as to take a value in the range [—1, 1]. When A0 ^ +1, then the local value 
of (pa is high whilst the value of (pb is low. Conversely, when A0 ^ —1, then the local (pa 
is low and (pb is high. The average order parameter values in Fig. 10 are: (a) (p = 0.25, 
(b) = 0.3, (c) = 1, (d) = 1.15, (e) = 1.2, (f) = 1.27. The most palpable 
change from the one component model is that the phase diagram is largely dominated by 
the striped profiles, with stripes appearing in the range 0.22 ^ ^ 1.28. Just outside the 
range of where bump structures are formed, we observe order parameter profiles which 



contain a mixture of bumps and stripes - see Fig. 10 a) - this value of must lie inside 



the coexistence region between the bump and stripe phases. As we increase the value of 
we enter the large region of parameter space where stripe structures are formed [Fig. [Io|b) 
and (c)], the only significant change as we increase from 0.3 to 1 is the decrease in the 
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width of the stripes; this is due to the fact that the typical length scale in the system is 



2t[ /km, where the wave number km given by Eq. (29), is inversely proportional to the average 
order parameter values 0^ and 0^. Increasing the value of further, we continue to observe 
striped profiles, but now there are points where the stripes of one species 'connect' to stripes 
of the other species - see Fig. [Io|(d) (these 'connections' appear as white lines in Fig. [lQ|(d)). 



Increasing further, we observe a mixture of holes and stripes [Fig. 10 e)]. Close to the 



instability curve Eq. (30) we find interesting profiles where we observe a mixture of stripes, 
holes and regions where the profile is approximately uniform 0„ 0^ [Fig. [lO|e)]. 
Various modulated structures are observed over a large range of parameter values. It would 
be possible to consider the structures formed for different values of the coupling coefficient 
7] and different values of the average order parameters, where 0a 7^ 0{,. However, here we do 
not make a systematic study of the entire parameter space and instead focus on the various 
bump formations. These structures closely resemble the configurations of particles/colloids 
in condensed matter systems and we believe that in this regime the model may be useful to 
understanding the fiuid and solid phases of such systems. 

B. Intermolecular interactions 

For the remainder of this paper, we pursue the idea that the bumps in this two component 
model represent two different types of molecules or colloidal particles suspended in a fiuid 
medium. We perform time simulations of the two component model choosing parameter 
values which result in the formation of bump structures. We run these simulations until 
the order parameter profiles reach an (almost) stationary state, which corresponds to being 
at (near) an energetic minimum. We then determine the coordinates of the particles by 
locating the position of the maximum of each of the peaks. The radial distribution functions 
are calculated by analysing these coordinates. We also calculate the effective pair potentials 



between the bumps. Later in Sec. |III C| we consider the nearest neighbour bond angles and 
the ordering in crystalline configurations. 

The bump phase in the two component model appears to behave in a similar manner to that 
of the one component model (cf. Fig. |6|. We observe bump structures when the average 
value of the order parameters 0a and 0{, are small. In particular, when (pa = 4>b = 4> and 
r = —0.9 we observe bumps within the range < < 0.15. We study and compare two 
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different systems: the symmetric case where qa = Qb = ^ and the asymmetric case where 
Qa = 1 and qb = 1.1. In the symmetric case, interactions between bumps of the same type 
{aa and bb) are identical in both components, but the nature of the interaction between a 
bump in (pa and a bump in 0^ (ab) is determined by the couphng term in the free energy. In 
the asymmetrical case, the different values of q mean that the size of the bumps are different 
in (pa and (pb, hence, all possible interactions aa, bb and ab are different. 
We begin by considering how the two order parameter profiles change as we alter their 
average values (pa = <Pb = <P- Thus we keep the concentration of the mixture fixed at c = 0.5, 
where 



(33) 



+ (pb 



We set the other parameter values to a = 1, r = —0.9 and 77 = 4. In Fig. 11 we display 



typical results. We plot the normalised difference between the two order parameters A(p 



(as defined in Eq. (32)). In the left column we show profiles from the symmetric case and 
in the middle column we display the profiles from the asymmetric system. In the right 
column we present the radial distribution functions, which are obtained by averaging over 
at least fifty runs, each with different realisations of the initial noise. The solid lines show 
the radial distribution functions for the symmetric case and the dashed lines show the 
asymmetric case. It is very apparent that this region of the parameter space shares many 
similarities with the one component model in both one and two dimensions. If we select a 
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small value of (p we find localised peaks surrounded by vacant shown in Fig. 

(a) and (b). We observe a tendency for bumps in (pa and in (pb to sit pairwise next to 
each other resembling configurations occurring in mixtures of oppositely charged colloidal 
particles [311436] . When (pa ~ (pb, the arrangement of the bumps also resembles snapshots 
of monovalent salts. It is very difficult to differentiate between the structures formed by 
the symmetric and asymmetric models for small values of (p. This is because structurally, 
there is very little difference between the two cases. If we examine the radial distribution 



functions gij (where i,j = a,b) for the symmetric and asymmetric systems [Fig. 11 'c)] we 
observe that the average distance between the different bumps seems to be independent of 
the q values (any differences between the curves is of the same order of magnitude as the 
statistical error). This is due to the large vacant areas, which means that there are not 
many bumps which are close to one another, especially between bumps in the same species 
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FIG. 11: In the left hand column we display typical simulation results for A(j) for the symmetrical case 
where Qa = qt = ^- In the middle column we display results from the asymmetrical case where = 1 
and qi, = 1.1. The orange regions show where there is a 4>a bump, while the blue show the 0^ bumps 
which are slightly smaller in the asymmetric mixture. In the right hand column the radial distribution 
functions gij (x) are shown for the symmetrical case (solid lines) and the asymmetrical case (dashed lines) . 
The parameter values are: a = 1, 77 = 4, r = —0.9, <j)a = 4)b = 4), where (a) - (c) ^ = 0, (d) - (f) ^ = 0.04, 
(g) - (i) = 0.06 and (j) - (1) = 0.15. 
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(aa and bb). 

As we increase the values of we find that the number of bumps of both species increases. 



In Fig. 11 d)-(f) we show the case where = 0.04 and in Fig. 11 g)-(i) we show the case 
where = 0.06. There is now a clear difference between the symmetric (d), (g) and the 
asymmetric (e), (h) cases. We observe a larger number of bumps in (pf, when qi, = 1.1. This 
is because the larger value of q reduces the length scale of the modulations, meaning that 
more bumps can be created before the value of 0f, becomes small (and negative) in vacant 
areas. There is an optimum value of (pa and (pf, in the vacant (uniform) areas which depends 
on the parameter values. This explains why increasing the value of (p increases the number 
of bumps (i.e., more modulations are needed in order to reach the optimum value of (p in the 
vacant regions). These intermediate values of (p produce profiles with bump configurations 
that resemble real fluid structures. However, in stark contrast to the one component system 
(shown in Fig. [7]^b)), we now find the formation of chains of alternating bumps reminiscent 
of structures observed in charged fluids. The radial distribution functions in Fig. [lT|f) and 
(i) show that the asymmetry induced by the different values of q begins to take effect at 
these intermediate values of (p. We observe that statistically the bumps sit closer together in 
the asymmetrical case, especially when two bumps in (pb are next to each other {bb, shown 
by green dashed line). This is due to the decreased size of the bumps in (ph, allowing them 
to sit slightly closer to their neighbours. 

Increasing the average order parameter values (p further we begin to observe the formation 
of crystalline structures as shown by Fig. [TT](j)-(l). The interesting thing is that now we 
observe square ordering of the particles instead of the hexagonal ordering which is present 
in the regular PFC model and the one component VPFC model. This implies that as we 
increase the concentration of one of the species from c = (almost a pure one-component 
system) to c = 0.5 there must be a transition from hexagonal to square ordering of bumps, 



this is something we return to below in Sec. Ill C Just as for the one component system, we 
find that there are more modulations in (pb when qb = 1.1. The profiles obtained with these 
parameter values resemble a compound crystal structure with vacancies and grain bound- 
aries. The radial distribution functions in Fig. [lT|^l) show that the smaller size difference 
of the (pb bumps in the asymmetric mixture has a large impact on the average position of 
the bumps in the structure compared to the symmetric mixture. This is because the higher 
concentration of particles forces them all closer together resulting in all pairs of bumps aa, 
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FIG. 12: Fits to the shape of individual bumps (cf. Fig. 8) in the (a) symmetric, when qi, = 1, and (b) 
asymmetric, when qb — 1.1, systems and the corresponding dips in the other order parameter profile which 
occur under the bumps in (c) the symmetric and (d) the asymmetric systems. In the symmetric profiles 
(a) and (c) the c/ja and (j)b curves are equal everywhere. The bump profile in (f>a is virtually identical in the 
cases where qi, — 1 and qi, = 1.1. These fits are then used to calculate the elTective pair potential between 
bumps, which are displayed in (e). The inset displays a magnification of the tails of Vij{x). The resulting 
pair potential V{x) between two bumps in (j)a when qi, — 1.1 lies on top of the aa,qb — 1 curve. The 
parameter values are: a — 1, r — —0.9, »7 = 4, and 4>a = 4'b = 0. 



hh and ah being closer together. 



In Fig. 12 a)-(b) we show the shape of the individual bumps in and obtained in the 



low density limit — )■ 0. To determine these radially symmetric profiles we fit functions of 



the form 9{x) as defined above in Eq. (21). The bumps in 0^ are virtually identical for both 
the symmetrical and asymmetrical systems. The (pa bump in the symmetric system and 
the 0f, bump in the asymmetric system decay to different values due to the different values 
of (f)a and (pii in the vacant areas of the asymmetrical system. We observe that in this two 
component model, a bump in one order parameter profile coincides with a small depression 
in the other order parameter profile. This is caused by the coupling term, which means that 
the combination of a bump in one order parameter and a hole in the other order parameter 
reduces the free energy of the system. In Fig. [I2|^c)-(d) we show the shape of the 'holes' 
which form in one profile under the bumps in the other order parameter field. These are 
determined the same way as the bump profiles: by fitting a function of the form shown in 



Eq. (21) to data points obtained from simulations. The depth of the holes is much smaller 
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in size than the height of the bumps. This is because the vacancy term prevents the hole 
from reaching large negative values of (pa oi (pb- 

Using the fitted functions shown in Fig. [I2|(a)-(d) we calculate effective pair potentials Vij (x) 
between the different particles in the system = a,b). We do this by determining the 
free energy for a system containing two bumps and their corresponding holes at various 
distances apart. In Fig. [I2|^e) we display the effective pair potentials for both the symmetric 
(solid lines) and the asymmetric (dashed lines) systems. The results show that there is 
an attraction between all of the bumps, just as we found for the one component system 
[Fig. |8](b)]. In both the symmetric and asymmetric cases we find that the attraction between 
two bumps from different species [ab) is stronger and occurs at a smaller value of x than 
that of two bumps of the same species {aa and bb). This explains the tendency for the 
bumps to form chains at intermediate values of 0a and (pi, [Figs. [ll1(d), (e), (g) and (h)] and 
square ordered crystalline structures at larger values of (pa and (pb [Fig. [TT|j)-(k)]. This is 
also consistent with the appearance of the large peak in Qabi^) which occurs at a smaller x 
value than the main peaks in gaa{x) and gbb{x) - see Figs. [TT|(c), (f ), (i) and (1). The effective 
pair potential Vaa{x) is almost identical in the symmetric and the asymmetric systems. This 
suggests that the small hole which appears in (pb has little effect on the interaction between 
the bumps. The major difference between the symmetrical and asymmetrical systems is 
that in the asymmetric mixture the minimum of the pair potentials Vab{x) and Vbb{x) are 
at smaller values of x than in the symmetric mixture. This is due to the reduced size of 
the (pb bumps in the latter. The minimum in Vbb{x) is at a slightly larger value of x than 
the minimum in Vaa{x) and the attraction is also much weaker (in fact it is so much weaker 
that the minimum is barely visible in this plot). This to some extent explains why the effect 
of the asymmetry is not visible for smaller values of (pa and 0^, but becomes apparent for 
larger values of (pa and 0^, where the vacant areas become smaller and we observe a close 
packing of the particles. 



C. Bond angles and the transition between hexagonal and square ordering 

In the two dimensional one component model (Eqs. ([T]) and (|4])) we observe hexagonally 
ordered structures for certain parameter values [Fig. |6](a) and Fig.[7]^c)]. However, in the two 
component model when (pa = <pb, we instead observe a square ordered crystalline structure 
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which alternates between species a and species b Figs.[ll[j)-(k)]. Thus, as the composition of 
the mixture is varied we should see a transition/ crossover from hexagonal to square ordering. 
The number of bumps observed in each field 0j depends on the respective average value 0,. 



When the concentration c ~ or c ~ 1, where c is defined in Eq. (33), (i.e., when either 
4'b ^ 4'a or (pa ^ (pb) then the resulting order parameter profile A0(x) has many more 
bumps of one type than of the other, and in these two limits we again observe hexagonal 



ordering. Note that c in Eq. (33) is not a bump concentration, but instead is a ratio between 
the two average order parameter values. As the 0j may take a negative value, for c = there 
are still a few bumps of a and similarly there are still some species b bumps when c = 1. 
When c = 0.5 the number of bumps is roughly the same in both species for the symmetrical 
case {qa = qb), but this is not necessarily true for the asymmetrical system (g^ 7^ qb). When 
4>a = 4>b and qa < qb there are more b bumps than a bumps. 



In Figs. 13 'a)-(c) we show the order parameter A0 for varying values of c for the symmetric 
mixture {qa = qb)- We fix the total 'density' (f)a+4>b = 0.24 and investigate how the crystalline 
structures change as the concentration c is varied. In Fig. [l3|^a), when c = we observe 
a profile which is dominated by species b bumps. The crystal is hexagonally ordered with 
some defects (these tend to occur in the vicinity of the (pa bumps). There are only a few 
(pa bumps, which means the bumps in b are usually sitting next to each other, resulting in 
them ordering themselves in a similar manner to that observed in the one component model 
Fig. [rj^c). Increasing the value of c from to 0.25, we observe a loss of crystalline structure, 
as shown in Fig. [l3Fb) . The loss of long range order is clearly visible in the associated radial 



distribution functions (not shown). The profile in Fig. 13 b) shows a somewhat amorphous 
structure which appears to include both square and hexagonal ordering in equal measure. 
Increasing the concentration further to c = 0.5, we observe a similar square ordering of 
bumps as in Figs. [II]^j)-(k). (in Fig. [ll|(j): qb = 1, whereas in Fig. [ll|(k): qb = 1.1, all other 



parameter values are the same). The crystalline structure in Fig. 13 c) at (pa = (pb = 0.12 



contains more vacancies and defects than the one in Fig. 11 j) at (pa = (pb = 0.15, as both 
average order parameter values are smaller. Owing to the symmetry induced by choosing 
qa = qb (i-e., 0a — ^ 0b as c — 7- 1 — c), a case with concentration c is equivalent to the case 
with concentration 1 — c. Thus Fig. [I3|^b) also shows the case c = 0.75 if one interchanges 
the orange and blue bumps. For this reason values c > 0.5 are not shown. 
For the asymmetric system the c — )■ 1 — c symmetry does not exist and we therefore show 
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FIG. 13: Plots of the order parameter Acj), in which bumps in species a appear in orange and bumps in 
species b appear in blue, for a constant total order parameter value (t>a + 4'b = 0.24. In (a) - (c) we show 
results from the symmetric mixture {qa = qt = where the concentration of species a is (a) c — 0, (b) 
c = 0.25 and (c) c = 0.5. In (d) - (h) we show the results from the asymmetric mixture {qa — 1 and 
qb = 1.1) where (d) c = 0, (e) c = 0.25, (f) c = 0.5, (g) c = 0.75 and (h) c = 1. The parameter values are: 
rj = A and r = —0.9. 



five cases for c varying from to 1 in Figs. 13 d) c = 0, (e) c = 0.25, (f) c = 0.5, (g) 
c = 0.75 and (h) c = 1. We again observe a transition from hexagonal ordering in Fig. [T3|^d) 
to square ordering in Fig. [I3|f) and back to hexagonal ordering in Fig. [I3|(h) as the value 
of c is increased from to 1. In between the highly structured states we observe the mixed 
ordered states [Figs. [l3|(e) and (g)] that were also present in the symmetrical system. By eye, 
it is very difficult to pick out the differences between the symmetrical and the asymmetrical 
cases. As previously discussed, the different value of in the asymmetrical system changes 
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FIG. 14: An example Delaunay triangulation is shown for a simple one-component case. In (a) we dis- 
play a typical order parameter profile for the one-component model where we observe isolated peaks. 
The coordinates of the maxima are calculated, these are shown as black points in (b). In (b) we show 
the Voronoi diagram (light blue polygon network) and the Delaunay triangulation (red triangles) for this 
particular set of coordinates. 

the shape, size and quantity of b bumps. In order to characterise and better understand the 
organisation of the crystalhne structures that are formed, we require a measure which may 
be used to quantify the structures and distinguish between hexagonal and square ordering in 
both the symmetric and the asymmetric systems. To do this, we use Delaunay triangulation 
[371 EH] to calculate the distribution of the bond angles p(9) between nearest neighbours. 
We could have used other measures from stochastic geometry |39j, which were used to 
characterise the hexagon- square transition in Benard convection |40j . 

The Delaunay triangulation is a triangulation of points (in our case the coordinates of the 
peaks of the bumps in both order parameter fields) which maximises the minimum angles of 
every triangle (i.e., avoids 'skinny' triangles). This triangulation can be calculated from the 
Voronoi diagram [HTJ EH] of any set of points on a 2d plane. The Voronoi diagram is a set of 
polygons, where each polygon represents an area in 2d space which is closer to a particular 
point than to any of the other points (i.e., the locus of points contained in each polygon 



is closer to the bump inside the polygon than any other bump). In Fig. 14 we show an 
example of how we calculate the Delaunay triangulation for a given order parameter profile. 
The example shows the triangulation for a one-component profile (as the pairing between 
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bumps in the two-component model makes the triangulation harder to see) but the process 
is apphed in the same manner to the two component modeL We take the coordinates of all 
the bumps to be our points on a 2d plane. We then calculate the Voronoi diagram (shown as 
the light blue lines in Fig. [l4|b)) which can be used to calculate the Delaunay triangulation 
(shown as the red lines in Fig. [I4|^b)). This can be done using any of the algorithms outlined 
in Refs. [371 EH]- For an efficient method of calculating Voronoi diagrams and Delaunay 
triangulations see Ref. [H]. (Note that Delaunay triangulation becomes degenerate when 
points appear in certain lines of symmetry. However, the initial noise added to the order 
parameter fields prevents bumps from forming in perfect symmetry). We use the statistics 
of the triangles in the Delaunay triangulation to characterise the structures produced by the 
bumps. 

We extract three quantities from the triangulation: the area of the triangles, the length of the 
sides and the angles in each of the triangles. This information is gathered for five different 
realisations of the initial noise profile for systems of size 200 x 200 and the information is 
sorted into bins. From these bins we obtain the probability distribution function for each 
quantity. Comparing the different distributions for various values of c allows us to observe 
how the triangles in the triangulation change as we go from hexagonal to square ordering. 
Here we concentrate on the probability distributions of the angles in the triangulation to 
characterise the crystalline structures. For results from the other measures see Ref. [12]. 



In Fig. 15 we display the probability distribution p{Q) for the triangle corner angles, as 
the concentration c is varied from to 1. We show results for the symmetric (solid red 
line) and the asymmetric (dashed blue line) systems. The distribution of the angles of the 
triangles clearly shows the transition between hexagonal and square ordering. When c = 
we observe hexagonal ordering in both the symmetric and the asymmetric systems, which 
results in the formation of roughly equilateral triangles in the Delaunay triangulation. This 
produces angle distributions which have a single peak at 60°, as shown in Fig. [iSj^a). As 
the value of c increases the structure changes to square ordering, transforming the triangles 
into right-angled triangles. Hence the angle distribution changes and we observe a peak 
slightly above the value 45° and another peak (half the size) slightly below 90°, as shown 
in Fig. [Tsj^c). Increasing the concentration further to c = 1 restores the hexagonal ordering. 



hence the angle distribution returns to the single peak at 60° [Fig. 15 e)]. In between the 
purely hexagonal and the purely square ordered structures we observe states where the 
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FIG. 15: The bond angle distribution p{Q) for a constant total order parameter value 4>a + <l>b = 0.24. 
The concentration of species a is (a) c = 0, (b) c = 0.25, (c) c — 0.5, (d) c — 0.75, and (e) c = 1 (cor- 



responding to the simulation snapshots shown in Fig. 13). Results for the symmetric system are shown 
as the red solid lines and the asymmetric system results are shown as blue dashed lines. The parameter 
values are: ga = 1, '/ = 4 and r = —0.9. 



distribution of bond angles is more evenly spread, with small peaks occurring just above 



45°, at around 60° and just below 90° [Figs. 15 b) and[l5[d)]. These represent the somewhat 
amorphous structures which lack the long range ordering which is present in the hexagonally 
and square ordered structures. The position of these peaks in the bond angle distributions 
p(0) does not depend on the quantity or size of the bumps and so the peaks occur in (almost) 
the same position for the symmetric and the asymmetric systems for all values of c (this is 
not the case for the area or length distributions). This makes the bond angle distributions 
ideal for comparing the structure of bump formations in different systems. On comparing 
the symmetrical and the asymmetrical cases we observe that p(0) appears smoother in the 
symmetrical case. The distribution function p(0) has a more jagged appearance for the 
asymmetrical mixture, which we believe is due to the fact that there are different sized 
bumps in this mixture, making it more difficult for the bumps to organise themselves into 



regular structures. In Figs. 15 a) and |15[e) the distributions appear very similar for the 
symmetric and the asymmetrical cases, however, in the other distributions (in particular, 
Figs. [Tsj^b) and[T5|^d)) we observe a distinct difference in the height of the three peaks. This 
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suggests that the transition between the different ordered states occurs differently in the 
symmetric and the asymmetric systems. 

To examine more closely the transition from the hexagonal to the square ordered states 
we introduce an order parameter $ which is calculated from the distribution of the angles 
from the Delaunay triangulation. We integrate the angle distributions over three regions 
which cover the three different peaks (these regions are determined arbitrarily from close 



examination of the angle distributions in Fig. 15 ) and define the quantities: 
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R2 = / p(e) dQ. (34) 
J72 

We then define the order parameter $ in the following way: 

* = (35) 

When a structure consists of mainly hexagonal configurations of bumps the value of $ is 
small (since $ — ?■ as -Rq and R2 — ?■ 0) and when a profile is dominated by square 
ordering the value of $ is large (since $ — )■ 00 as -Ri — )■ 0). Calculating this quantity for the 
angle distributions for different values of c gives us a measure for the hexagonal vs. square 
ordering of the bumps. 

In Fig. [16] we show how the order parameter $ changes with the concentration c for the 
symmetric (solid red line) and the asymmetric (dashed blue line) mixtures. Both curves show 
a smooth continuous transition from hexagonal ordering to square ordering and back again. 
The different sized bumps in the asymmetric system break the symmetry around c = 0.5 
and we observe that the maximum (which corresponds to the strongest square ordering) 
occurs at around c ~ 0.6, and is actually higher than the peak in the symmetric case. The 
transition to and from square ordering appears to be slightly sharper in the asymmetrical 
case. Even though there is a difference in the transition between the different ordered states 
in the symmetric and the asymmetric mixtures, they appear to be qualitatively similar. It 
may be the case that for a larger difference in the values of and qb a different type of 
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FIG. 16: Plot showing the order parameter $ (defined in Eq. (35)) as a function of the concentration 
c of species a, where (pa + (f'b = 0.24. The red solid line and points show the symmetric case and the 
blue dashed line and points show the asymmetrical case. The parameter values are: qa = 1, rj ^ A and 
r = -0.9. 



transition from hexagonal to square ordering might occur, e.g., a discontinuous transition. 
However, the effect of varying the ratio Qa/qb is not studied in detail here. 



IV. CONCLUSIONS 



In this paper we have investigated the VPFC model and its application to materials mod- 
elling. We first considered the one-component model proposed by Chan et al. in Ref. [T3] . 
We determined the linear stability of the homogeneous state and discussed the dispersion 
relation. We examined the phase behaviour in one dimension and calculated the phase dia- 
gram, computing exactly the location of the tricritical point between the homogeneous and 
periodic states, and identified the region of phase space where localised structures occur. 
Focusing on the latter region of the phase diagram, we investigated the localised steady 
state profiles and discussed the slanted homoclinic snaking which occurs in the bifurcation 
diagrams. The one-component model was also studied in two dimensions and we deter- 
mined the phase diagram, radial distribution functions and effective pair potentials from 
our simulation data. Some of the behaviour we have identified - the presence of transitions 
resembling transitions from a solid phase to a liquid phase and then to a gas-like phase 
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- replicates behaviour observed in nonconserved systems [13]. In section III of the paper, 
we extended the model to include two coupled order parameter fields. We have considered 
how the coupling affects the linear stability of flat fllms and then briefly touched on the 
phase behaviour of this two component model. We have focused on the bump structures 
which form, considering both a symmetrical mixture where the bumps are of equal size and 
an asymmetrical system where one of the bump species is slightly smaller than the other 
species. The radial distribution functions and effective pair potentials for these systems are 
somewhat similar to those in binary mixtures of oppositely charged colloidal particles. We 
have investigated how varying the concentration c of the mixture produces a crossover from 
hexagonal to square ordered crystalline structures and how the transition differs between 
the symmetrical and the asymmetrical systems. 

A key issue on which we should comment concerns the question of what precisely does the 
order parameter proflle 0(x, t) in the VPFC model represent? In the regular PFC model, 
the phase with the uniform flat proflle is taken to represent the liquid phase, whilst the 
bump phase corresponds to the crystalline solid. This interpretation is underpinned by the 
fact that the regular PFC can be derived from density functional theory (DFT) [2] and 
dynamical density functional theory (DDFT) [10], which is a theory for the dynamics of 
a system of interacting Brownian (colloidal) particles [lH-STj. DFT [l8] - [50] is a statistical 
mechanics theory for the one-body number density p(x) of a system of particles, where 
p(x) = (p(x)) and where p(x) = ~ ^i) density operator and (■) denotes a 

statistical ensemble average [IB] . The central quantity in DFT is the Helmholtz free energy 
functional F[p\ and the equilibrium fluid density proflle p(x) is that which minimises the 
grand free energy Q[p] = F[p]— ^ J dxp(x). The DDFT for Brownian particles [HHH] takes 
as input this functional and so yields the correct equilibrium fluid density proflle. Making 
a truncated gradient expansion approximation for F[p], expanding the free energy around 
that of a reference liquid state with uniform density po, one can argue that the free energy is 
approximately given by Eqs. ^ and ([s]), where the order parameter 0(x) oc p(x) — po- Thus 
it is clear that in a bulk liquid, where p(x) is a constant, so too is 0(x) a constant and in the 
solid phase, where p(x) consists of a periodic array of density peaks, then 0(x) also contains 
periodic modulations. However, there are some problems extending this interpretation to the 
VPFC. Consider for example Fig. [7](a) where we see a few isolated localised peaks surrounded 
by a uniform background where 0(x) ^ 0. Maintaining the above PFC interpretation, this 
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would correspond to a few individual 'frozen' particles, surrounded by a fluid of mobile 
particles. One might be tempted to think of this as some sort of glass transition [SUHS2], 
but the glass transition is a collective phenomenon: in a glass one does see 'dynamical 
heterogeneity' i.e., regions where the particles are totally jammed and other regions which 
are more mobile, but to our knowledge one never sees a single particle that is jammed on 
its own surrounded by more mobile particles. Thus, it may be possible to assume this 
interpretation may be maintained for the VPFC, i.e., by considering the localised peaks 
surrounded by a uniform background to be a dynamically heterogeneous glassy system, but 
there are problems with this point of view. 

An alternative interpretation for the order parameter profile in the VPFC model is that 
4>{x) is related to a coarse-grained density profile (rather than an ensemble average density 
profile) for the system p(x, t), i.e., 0(x, t) oc p(x, t). Following Ref. |17], we may define 
the temporally coarse-grained density profile for a system of Brownian colloidal particles 
as p(x, t) = J K{t — t')p{'x,t)dt' , where K{t) is a normalised function of finite support 
which defines a time window over which the density is coarse-grained. One can then argue 
[17j that the time evolution equations for p(x, t) must be very similar or even the same 
as the DDFT equations for the time evolution of the ensemble average density p(x, t), as 
long as the width in time r of K{t) is large enough. By choosing the time r so that it is 
large compared to the time between the colloidal particles receiving Brownian 'kicks' from 
the solvent, but is short compared to the diffusive time scale, corresponding to the typical 
time for a particle to diffuse a distance equal to its own diameter, then the coarse-grained 
density p(x, t) and the order parameter (/)(x, t) will be quantities which contain peaks, each 
of which correspond to an individual particle in the system. Thus, in a low density colloidal 
suspension one should see isolated peaks in the coarse-grained density, surrounded by regions 
where 0(x, t) ^ 0, corresponding to no particles being present in that region of the system. 
This is the justification for the interpretation made by Chan et al. in Ref. [Uj, that the peaks 
in the order parameter correspond to particles and the uniform background corresponds to 
a portion of solvent free of particles. In order to observe the long time Brownian motion of 
the particles in this description, one should add a stochastic noise term to the dynamical 



equations for the system (23), that continuously drives the system (as opposed to the small 
amount of noise that is present in our initial order parameter profiles). However, in numerical 
simulations there can be problems with such an approach, because the particles can become 
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pinned in place by the discrete grid on which they are defined, and so do not move. We did 
not make a detailed investigation of the of the VPFC model with additional noise. Further 
issues arise as the noise renormalises the parameters of the continuum model. 
There are state points in the PFC and VPFC phase diagram where all possible interpreta- 
tions of (j) break down: these are the state points where the equilibrium state is the stripe 



or the hole phase, such as those displayed in Fig. [TOj Systems of spherical particles do 
not have an ensemble average density profile p nor a coarse-grained density profile p with 
stripes/holes, unless the particles in the system interact via pair potentials containing com- 
peting attractive and repulsive parts [TSl ESI [51] • We must conclude that for the parameter 
values corresponding to these state points, the gradient expansion that is implicit in the 
PFC and VPFC free energy functionals has broken down and that these order parameter 
profiles are unphysical. 

The radial distribution functions for the one-component model displayed in Fig. [T] (see also 
Fig. 5 of Ref. [13]) are very similar to those in real fiuids. We observe static correlations 
which are very similar to what one observes in fiuids. Increasing the value of increases the 
number of bumps and close packing causes long range (crystalline) ordering of the bumps. 
Calculating the effective pair potential between isolated pairs of bumps, we find a pair 
potential having an attractive minimum at a pair separation distance which is slightly larger 
than the diameter of the bumps. Thus, the interactions and correlations between bumps 
share certain features with some colloidal fiuids |H3] . We also extend the model to consider a 
two component mixture, with a simple repulsive coupling between the two order parameter 
profiles. At low values of the bumps commonly appear in pairs and at intermediate values 
they tend to form chains. At higher values of the system exhibits crystalline ordering. The 
appearance of these structures is somewhat reminiscent of the arrangement of the particles 
in a binary mixture of oppositely charged colloidal particles - see e.g., Ref. [Ml [56] and 
references therein. The radial distribution functions and the effective pair potentials show 
there is a fairly strong attraction between bumps of the opposite species a and h. The 
minimum in the ah effective pair potential is at a shorter pair separation distance than 
the minimum in the aa and hh pair potentials and so we observe square ordering when 
the concentration c ~ 1/2 and is high enough for the bumps to pack into a crystalline 
structure. However, when c ~ or c ~ 1, we observe hexagonal ordering and so we observe 
a transition from hexagonal to square ordering as the concentration c is varied. We find 
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that this transition occurs smoothly but can become skewed by changing the size of one of 
the species of bumps (g^ 7^ gt). 

It would be interesting to further investigate the effect that varying the ratio qal^b has on 
the structures which form. In particular, determining the range of values of qa and for 
which bump profiles form in the 2d system would allow one to determine the range of size 
ratios of particles (bumps) that can be modelled. The transition between hexagonal and 
square structures could then be studied for systems with very different sized bumps and 
if the VPFC in this regime continues to be able to model mixtures of charged colloidal 
particles, then a wide range of different crystal structures should be observed [31]. 
Note also that the localised structures that we observe are not a unique property of the 
VPFC model but are in fact also present in the regular PFC model for a small range of 
parameter values outside the limit of linear stability. This is something that we will focus 
on in future work. 
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